%MACRO cox_analysis(in);

*****************************************************************************
* Program: COX_ANALYSIS					   						
* 																			
* Description: Perform Cox reegression and related plots including
*				and proportionality checks.
*						
* Last modified: September 6, 2015												
* Created by: Martin Wegman													
*****************************************************************************;

%let din = cox1; * From propensity.sas;

data loc.&din;
	set loc.&in;
run;

options PAPERSIZE=LETTER
	orientation=portrait
	topmargin="0.00in"
	bottommargin="0.00in"
	leftmargin="0.00in"
	rightmargin="0.00in"
	nodate nonumber;

options nodate nonumber;
ods noproctitle;
ODS ESCAPECHAR='^'; ;
%let date = %sysfunc(date(),DATE7.);

****************************************************************************;
****************************************************************************;
************************     Model variables       *************************;
****************************************************************************;
****************************************************************************;

* Vars after treatment assignment;

* Continuous vars;
%let var3 =  cs01_craving so_rec_d so_amb_d so_ts_d; 
* Categorical vars;
%let var4 = ;
****************************************************************************;
****************************************************************************;
*****************             Parent macros           **********************;
****************************************************************************;
****************************************************************************;

proc format library = loc;
value combin2_ 	1="CDDC(SR) + VTC(SR)"
				2="CDDC(SV14)+ VTC(SR)"
				3="CDDC(SV0) + VTC(SR)"
				4="CDDC(SR) + VTC(SV90)"
				5="CDDC(SV14) + VTC(SV90)"
				6="CDDC(SV0) + VTC(SV90)"
				7="CDDC(SR) + VTC(SV30)"
				8="CDDC(SV14) + VTC(SV30)"
				9="CDDC(SV0) + VTC(SV30)";
run;

%MACRO weighted_curves(var,type,site,data,title1,title2,lab,plot_type);

ods exclude all;
proc sort data = loc.&data; 
	by combination _imputation_;
proc phreg data=loc.&data;
	by combination _imputation_; 
	model T&type._&var*Y&type._&var(0)=&var3; 
	strata &site; 
	freq wt2/notruncate; 
	baseline out=temp survival=S LOGLOGS=LS2;
run;
ods exclude none;

proc sort data = temp; 
	by &site combination T&type._&var;
PROC MEANS DATA=temp noprint;
	by &site combination;
  	CLASS T&type._&var;
  	VAR &plot_type;
  	OUTPUT OUT=temp2 MEAN= ;
RUN;

data loc.temp3;
	set temp2;
	if _TYPE_ = 1;
	drop _TYPE_ _FREQ_;
run;

%if &plot_type = LS2 %then %do;

* Export for work in R;
libname out xport "&lib\mydataY&type._&var..xpt";
data out.mydata;
	set loc.temp3;
	combin =combination;
	drop combination;
run;

%end;

title &title1;
title2 "&title2";

proc sgpanel data=loc.temp3;
	panelby combination /COLUMNS= 3 NOVARNAME;
	*refline 50 / AXIS= X;
	scatter y= &plot_type x= T&type._&var /group=&site;
	LOESS y= &plot_type x= T&type._&var /group=&site DEGREE=2;
	ROWAXIS label = "&lab";
	COLAXIS label ="Time from release/discharge (days)";
	format combination combin2_.;
	label &site = "Group";
run; 
quit;

title;
title2;

%MEND;

%MACRO master_weighted_curves(var,type_list,site,data,plot_type);

%local i;
%let i = 1;
%let vals = a b c d e f;

%if &plot_type = LS2 %then %do;
	%if &var =opi %then %let text  = 'Figure S6$. Log of weighted, adjusted cumulative hazard of opiate relapse.';
	%else %if &var =amp %then %let text  =  'Figure S7$. Log of weighted, adjusted cumulative hazard of amphetamine relapse.';
	%else %if &var =bzd %then %let text = 'Figure S. Log of weighted, adjusted cumulative hazard of benzodiazepine relapse.';
	%else %if &var =any %then %let text  = 'Figure S8$. Log of weighted, adjusted cumulative hazard of illicit-drug relapse.';
	%else %if &var =mtd %then %let text  = 'Figure S9$. Log of weighted, adjusted cumulative hazard of methadone linkage.';
	%else %if &var =bup %then %let text  = 'Figure S. Log of weighted, adjusted cumulative hazard of buprenorphine use.';
	%else %if &var =hiv %then %let text  =  'Figure S. Log of weighted, adjusted cumulative hazard of HIV infection.';
%end;
%else %do;
	%if &var =opi %then %let text  = 'Figure S14$. Weighted, adjusted survival function for opiate relapse.';
	%else %if &var =amp %then %let text  =  'Figure S15$. Weighted, adjusted survival function for amphetamine relapse.';
	%else %if &var =bzd %then %let text = 'Figure S. Weighted, adjusted survival function for benzodiazepine relapse.';
	%else %if &var =any %then %let text  = 'Figure S16$. Weighted, adjusted survival function for illicit-drug relapse.';
	%else %if &var =mtd %then %let text  = 'Figure S17$. Weighted, adjusted survival function for methadone linkage.';
	%else %if &var =bup %then %let text  = 'Figure S. Weighted, adjusted survival function for buprenorphine use.';
	%else %if &var =hiv %then %let text  =  'Figure S. Weighted, adjusted survival function for HIV infection.';
%end;

%if &plot_type = LS2 %then %do;
	%if &var = opi %then %let lab  = Complementary log-log of the probability of no opiate use;
	%else %if &var =amp %then %let lab  = Complementary log-log of the probability of no amphetamine use;
	%else %if &var =bzd %then %let lab = Complementary log-log of the probability of no benzodiazepine use;
	%else %if &var =any %then %let lab  = Complementary log-log of the probability of no illicit-drug use;
	%else %if &var =mtd %then %let lab  = Complementary log-log of the probability of no methadone linkage;
	%else %if &var =bup %then %let lab  = Complementary log-log of the probability of no buprenorphine use;
	%else %if &var =hiv %then %let lab  =  Complementary log-log of the probability of no HIV infection;
%end;
%else %do;
	%if &var = opi %then %let lab  = Probability of no opiate use;
	%else %if &var =amp %then %let lab  = Probability of no amphetamine use;
	%else %if &var =bzd %then %let lab = Probability of no benzodiazepine use;
	%else %if &var =any %then %let lab  = Probability of no illicit-drug use;
	%else %if &var =mtd %then %let lab  = Probability of no methadone linkage;
	%else %if &var =bup %then %let lab  = Probability of no buprenorphine use;
	%else %if &var =hiv %then %let lab  =  Probability of no HIV infection;
%end;

%do %while("%scan(&type_list,&i)" NE "");

	%let type = %scan(&type_list, &i);
	%let val = %scan(&vals, &i);
	%let title1 = %sysfunc(translate(&text,&val,$));

	%if &type =A %then %let text2 = (Missing measurements are ignored);
	%else %if &type =B %then %let text2 = (Missing measurements coded as an event);
	%else %if &type =C %then %let text2 = (Missing measurements coded as a non-event);
	%else %if &type =D %then %let text2 = D- Worst case scenario;
	%else %if &type =E %then %let text2 = E- Missing equivalent to censored;
	%else %if &type =F %then %let text2 = F- Intervening missing a non-event/trailing missing an event;
	%else %if &type =G %then %let text2 = (Reverse event coding where missing is ignored);
	%else %if &type =H %then %let text2 = (Missing measurements coded with differential risk of event);

	%weighted_curves(&var,&type,&site,&data,&title1,&text2,&lab,&plot_type);

	%let i = %eval(&i + 1);

%end;

%if &plot_type = LS2 %then %do;
%if "&var" ne "mtd"  %then %do;

%let text = 
Note: For each stage 1 imputation combination, 
plots of the log of the cumulative hazard of the event 
estimated using Cox regression models stratified by group, 
weighted by the inverse propensity of seeking treatment, 
and adjusted for measurements of opiate cravings 
and ambivalence, recognition, and taking steps towards change. 
Estimates are simple averages of stage 2 imputation results (M=5).
Local regression smoothing lines are also provided for ease of interpretation. 
The heading for each subplot reads CDDC imputation model + 
VTC imputation model.
Each set of plots uses a separate event definition as described in the subtitle.^2n;

ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}&text";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}VTC: Voluntary Treatment Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}CDDC: Compulsory Drug Detention Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SR: Stochastic regression imputation";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV90: Single value imputation of length of inpatient stay (90 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV30: Single value imputation of length of inpatient stay (30 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV14: Single value imputation of time of first outcome measurement (14 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV0: Single value imputation of time of first outcome measurement (0 days)";

%end;
%end;

%else %do;
%if "&var" ne "mtd"  %then %do;

%let text = 
Note: For each stage 1 imputation combination, 
plots of the survival function 
estimated using Cox regression models stratified by group, 
weighted by the inverse propensity of seeking treatment, 
and adjusted for measurements of opiate cravings 
and ambivalence, recognition, and taking steps towards change. 
Estimates are simple averages of stage 2 imputation results (M=5).
Local regression smoothing lines are also provided for ease of interpretation. 
The heading for each subplot reads CDDC imputation model + 
VTC imputation model.
Each set of plots uses a separate event definition as described in the subtitle.^2n;

ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}&text";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}VTC: Voluntary Treatment Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}CDDC: Compulsory Drug Detention Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SR: Stochastic regression imputation";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV90: Single value imputation of length of inpatient stay (90 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV30: Single value imputation of length of inpatient stay (30 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV14: Single value imputation of time of first outcome measurement (14 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV0: Single value imputation of time of first outcome measurement (0 days)";
%end;
%end;



%MEND;


****************************************************************************;
****************************************************************************;
*  Plot weighted cum hazard curves to assess proportionality assumption   **;
****************************************************************************;
****************************************************************************;

%let file = "&out_dir.\CDDC Urine Supp Cox CLL CumHaz &date..rtf";
ods rtf file = &file Style = custom_narrow startpage = Yes BODYTITLE;

%let plot_type = LS2;
%let data = cox1;
%let site = site3;
%let var = opi;
%let type_list = A B C H;

%master_weighted_curves(&var,&type_list,&site,&data,&plot_type);

%let plot_type = LS2;
%let data = cox1;
%let site = site3;
%let var = amp;
%let type_list = A B C H;

%master_weighted_curves(&var,&type_list,&site,&data,&plot_type);

%let plot_type = LS2;
%let data = cox1;
%let site = site3;
%let var = any;
%let type_list = A B C H;

%master_weighted_curves(&var,&type_list,&site,&data,&plot_type);

%let plot_type = LS2;
%let data = cox1;
%let site = site3;
%let var = mtd;
%let type_list = A B C;

%master_weighted_curves(&var,&type_list,&site,&data,&plot_type);

%let plot_type = LS2;
%let data = cox1;
%let var = mtd;
%let type = G;
%let site = site3;
%let title1 = %str(Figure S10d. Log of weighted, adjusted cumulative hazard of no methadone-retention.);
%let title2 = (Reverse event coding where missing is ignored);
%let lab = Complementary log-log of the probability of methadone retention;

%weighted_curves(&var,&type,&site,&data,&title1,&title2,&lab,&plot_type);

%let text = 
Note: For each stage 1 imputation combination, 
plots of the log of the cumulative hazard of the event 
estimated using Cox regression models stratified by group, 
weighted by the inverse propensity of seeking treatment, 
and adjusted for measurements of opiate cravings 
and ambivalence, recognition, and taking steps towards change. 
Estimates are simple averages of stage 2 imputation results (M=5). 
For methadone-retention time, the event was defined as
the first urine measurement for which the participant 
did not test positive for methadone.
Local regression smoothing lines are also provided for ease of interpretation. 
The heading for each subplot reads CDDC imputation model + 
VTC imputation model.
Each set of plots uses a separate event definition as described in the subtitle.^2n;

ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}&text";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}VTC: Voluntary Treatment Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}CDDC: Compulsory Drug Detention Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SR: Stochastic regression imputation";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV90: Single value imputation of length of inpatient stay (90 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV30: Single value imputation of length of inpatient stay (30 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV14: Single value imputation of time of first outcome measurement (14 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV0: Single value imputation of time of first outcome measurement (0 days)";

ods rtf close;


****************************************************************************;
****************************************************************************;
***********   Figure 3 - Adjusted Survival Curves for Opiates    ***********;
****************************************************************************;
****************************************************************************;

%let site = site3;
%let data = cox1;

ods exclude all;
proc sort data = loc.&data; 
	by combination _imputation_;
proc phreg data=loc.&data;
	where combination = 1;
	by combination _imputation_; 
	model TA_opi*YA_opi(0)=&var3; 
	strata &site; 
	freq wt2/notruncate; 
	baseline out=temp survival=S LOGLOGS=LS2;
run;
ods exclude none;

proc sort data = temp; 
	by &site combination TA_opi;
PROC MEANS DATA=temp noprint;
	by &site combination;
  	CLASS TA_opi;
  	VAR S;
  	OUTPUT OUT=temp2 MEAN= ;
RUN;

data loc.temp3;
	set temp2;
	if _TYPE_ = 1;
	if &site = 0 then group = 1;
	else if &site = 1 then group = 0;
	drop _TYPE_ _FREQ_;
run;

proc sort data = loc.temp3; by group; run;

options papersize=("8in", "6in")
	orientation=landscape
	topmargin=".5in"
	bottommargin=".25in"
	leftmargin=".25in"
	rightmargin=".25in"
	nodate nonumber;

%let var = opi;
%let type = A;
%let file = "&out_dir.\CDDC Urine Figure 3 &date..pdf";
ods noproctitle escapechar='^'; 
ods pdf file = &file style= custom_narrow_pdf2;



proc sgplot data=loc.temp3;
	step y= S x= TA_opi /group=group ;
	label S = "Probability of no opiate use";
	label TA_opi ="Time from release/discharge (days)";
	format group site2_.;
	XAXIS VALUES= (0 30 90 180 270 365);
	title "Figure 3. Adjusted probability of no opiate use, by group.";
	KEYLEGEND /TITLE= "Group" ACROSS= 1 LOCATION = INSIDE POSITION= TOPRIGHT;
run; 
quit;

title;

ods pdf close;

options PAPERSIZE=LETTER
	orientation=portrait
	topmargin="0.00in"
	bottommargin="0.00in"
	leftmargin="0.00in"
	rightmargin="0.00in"
	nodate nonumber;

****************************************************************************;
****************************************************************************;
*****************      Adjusted Survival Curves       **********************;
****************************************************************************;
****************************************************************************;

options nodate nonumber;
%let file = "&out_dir.\CDDC Urine Supp Adj survival curves &date..rtf";
ods rtf file = &file Style = custom_narrow startpage = Yes BODYTITLE;

%let plot_type = S;
%let data = cox1;
%let site = site3;
%let var = opi;
%let type_list = A B C H;

%master_weighted_curves(&var,&type_list,&site,&data,&plot_type);

%let plot_type = S;
%let data = cox1;
%let site = site3;
%let var = amp;
%let type_list = A B C H;

%master_weighted_curves(&var,&type_list,&site,&data,&plot_type);

%let plot_type = S;
%let data = cox1;
%let site = site3;
%let var = any;
%let type_list = A B C H;

%master_weighted_curves(&var,&type_list,&site,&data,&plot_type);

%let plot_type = S;
%let data = cox1;
%let site = site3;
%let var = mtd;
%let type_list = A B C;

%master_weighted_curves(&var,&type_list,&site,&data,&plot_type);

%let plot_type = S;
%let data = cox1;
%let var = mtd;
%let type = G;
%let site = site3;
%let title1 = %str(Figure S18d. Weighted, adjusted survival function for methadone-retention.);
%let title2 = (Reverse event coding where missing is ignored);
%let lab = Probability of methadone retention;

%weighted_curves(&var,&type,&site,&data,&title1,&title2,&lab,&plot_type);

%let text = 
Note: For each stage 1 imputation combination, 
plots of the survival function 
estimated using Cox regression models stratified by group, 
weighted by the inverse propensity of seeking treatment, 
and adjusted for measurements of opiate cravings 
and ambivalence, recognition, and taking steps towards change. 
Estimates are simple averages of stage 2 imputation results (M=5).
For methadone-retention time, the event was defined as
the first urine measurement for which the participant 
did not test positive for methadone.
Local regression smoothing lines are also provided for ease of interpretation. 
The heading for each subplot reads CDDC imputation model + 
VTC imputation model.
Each set of plots uses a separate event definition as described in the subtitle.^2n;

ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}&text";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}VTC: Voluntary Treatment Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}CDDC: Compulsory Drug Detention Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SR: Stochastic regression imputation";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV90: Single value imputation of length of inpatient stay (90 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV30: Single value imputation of length of inpatient stay (30 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV14: Single value imputation of time of first outcome measurement (14 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV0: Single value imputation of time of first outcome measurement (0 days)";

ods rtf close;


****************************************************************************;
****************************************************************************;
*****************        Print Hazard ratios          **********************;
****************************************************************************;
****************************************************************************;
proc format;

value group_ 	1 = "0 days < Time < 50 days"
				2 = "Time = 90 days"
				3 = "Time = 180 days"
				4 = "Time = 270 days"
				5 = "Time = 365 days";

 picture lc (round) /* Lower limit requires a left parenthesis and a comma on the right */
 			low-high = '009.99,'
 			( prefix = '(' );
 picture uc (round)
 			low-high = '009.99)'; /* Upper limit requires only a right parenthesis */
run;

%MACRO HR_estimates(data,var,type,site,ref,cox_type,file);

%let T = T&type._&var;
%let Y = Y&type._&var;

ods exclude all;
%if &cox_type = 1 %then %do;

	proc phreg data = loc.&data;
		by combination _IMPUTATION_;
		class &site (ref=&ref);
		format &site site3b_.;
		model &T*&Y(0)= &site /ties=efron;
		ods output  ContrastEstimate=est;
		contrast "Site" &site 1 /ESTIMATE=Parm;
	run;

%end;
%else %if &cox_type = 2 %then %do;

	proc phreg data = loc.&data;
		by combination _IMPUTATION_;
		class &site (ref=&ref);
		format &site site3b_.;
		weight wt2;
		model &T*&Y(0)= &site /ties=efron;
		ods output  ContrastEstimate=est;
		contrast "Site" &site 1 /ESTIMATE=Parm;
	run;

%end;
%else %if &cox_type = 3 %then %do;

	proc phreg data = loc.&data;
		by combination _IMPUTATION_;
		class &site (ref=&ref);
		format &site site3b_.;
		weight wt2;
		model &T*&Y(0)= &site &var3 /ties=efron;
		ods output  ContrastEstimate=est;
		contrast "Site" &site 1 /ESTIMATE=Parm;
	run;

%end;
%else %if &cox_type = 4 %then %do;
	
	proc phreg data = loc.&data;
		by combination _IMPUTATION_;
		class &site (ref=&ref);
		format &site site3b_.;
		weight wt2;
		model &T*&Y(0)= &site site50lt &var3 /ties=efron;
		if &T > 50 then site50lt = &site*log((&T-50));
		else site50lt = 0;
		ods output  ContrastEstimate=est;
		contrast "Site" &site 1 /ESTIMATE=Parm;
		contrast "Sitea" &site 1 site50lt 3.688 /ESTIMATE=Parm;
		contrast "Siteb" &site 1 site50lt 4.868 /ESTIMATE=Parm;
		contrast "Sitec" &site 1 site50lt 5.394 /ESTIMATE=Parm;
		contrast "Sited" &site 1 site50lt 5.753 /ESTIMATE=Parm;
	run;

%end;


proc sort data=est;
	by combination contrast _imputation_;
run;

proc mianalyze data=est;
	title;
	by combination contrast;
	modeleffects estimate;
	stderr stderr;
	ods output ParameterEstimates=parms;
run;
ods exclude none;

data &var._&type&cox_type;
	set parms;

	if Parm = "estimate";

	HR = exp(-Estimate);
	LCL = Estimate - 1.96*StdErr;
	UCL = Estimate + 1.96*StdErr;
	lowerCI = exp(-UCL);
	upperCI = exp(-LCL);
	CI = "(" || put(lowerCI,z5.3) || " - " || put(upperCI,z5.3) || ")";

	if contrast = "Site" then group = 1;
	else if contrast = "Sitea" then group = 2;
	else if contrast = "Siteb" then group = 3;
	else if contrast = "Sitec" then group = 4;
	else if contrast = "Sited" then group = 5;

	keep group combination HR CI LCL UCL LCLMean UCLMean lowerCI upperCI;
run;

%MEND;

%MACRO master_HR_estimates(var,type_list,data,site,ref);

	%local i j;
	%let i = 1;
	%let cox_type_list = 1 3 4;

	*Compute HRs;

	%do %while("%scan(&type_list,&i)" NE "");
		%let j = 1;
		%do %while("%scan(&cox_type_list,&j)" NE "");
			%let type =	%scan(&type_list,&i);
			%let cox_type = %scan(&cox_type_list,&j);
			%HR_estimates(&data,&var,&type,&site,&ref,&cox_type,&file);
			%let j = %eval(&j + 1);
		%end;

		%let i = %eval(&i + 1);
	%end;

%MEND;

%MACRO print_HR_estimates(var,type_list,title1);

	* Print HRs;
	%local i j;
	%let i = 1;
	%let cox_type_list = 1 3 4;
	%let vals = a b c d e f;
	%let vals2 = i ii iii iv;

	%do %while("%scan(&type_list,&i)" NE "");

		%let val = %scan(&vals, &i);
		%let type =	%scan(&type_list,&i);
		%let j = 1;

		%do %while("%scan(&cox_type_list,&j)" NE "");

			%let cox_type = %scan(&cox_type_list,&j);
			%let val2 = %scan(&vals2, &j);

			%if &cox_type =1 %then %let text = Unadjusted;
			%else %if &cox_type =3 %then %let text = Adjusted;
			%else %if &cox_type =4 %then %let text = %nrstr(Adjusted, time-varying);
			
			%let text1 = %sysfunc(TRANWRD(&title1,$,&val));
			%let text1 = %sysfunc(TRANWRD(&text1,#,&val2));
			%let text1 = %sysfunc(TRANWRD(&text1,@,&text));

			%if &type =A %then %let text2 = (Missing measurements are ignored);
			%else %if &type =B %then %let text2 = (Missing measurements coded as an event);
			%else %if &type =C %then %let text2 = (Missing measurements coded as a non-event);
			%else %if &type =D %then %let text2 = D- Worst case scenario;
			%else %if &type =E %then %let text2 = E- Missing equivalent to censored;
			%else %if &type =F %then %let text2 = F- Intervening missing a non-event/trailing missing an event;
			%else %if &type =G %then %let text2 = (Reverse event coding where missing is ignored);
			%else %if &type =H %then %let text2 = (Missing measurements coded with differential risk of event);

			%if (&cox_type = 1
				OR &cox_type = 2
				OR &cox_type = 3) %then %do;
				title "&text1";
				title2 "&text2";
				proc print data = &var._&type&cox_type label noobs; 
					var Combination /style(header data)={just=l};
					var HR CI /style(header data)={just=c};
					label combination="Combination";
					label HR="Hazard Ratio";
					label CI="95% Confidence Interval";
					format HR z5.3;
				run;
				title; title2;
			%end;
			%else %if &cox_type = 4 %then %do;
				title "&text1";
				title2 "&text2";
				proc tabulate data= &var._&type&cox_type;
					class combination group;
					var HR lowerCI upperCI;
					table 	Combination,
							group=' ' *(HR*max='' * [style=[just=C cellwidth=75] f=z5.3]
 									lowerCI*max="" * [style=[just=R cellwidth=65 ] f=lc.]
 									upperCI*max="" * [style=[just=L cellwidth=65 ] f=uc.]);	
					format group group_.;
					label combination="Combination";
					label HR="Hazard Ratio";
					label lowerCI ='95%';
					label upperCI ="C.I.";
				run;
				title; title2;
			%end;

			
			%let j = %eval(&j + 1);

		%end;

		%let i = %eval(&i + 1);

	%end;

%if "&var" ne "mtd"  %then %do;

%let text = 
Note: For each stage 1 imputation combination, 
hazard ratios and 95% confidence intervals 
comparing the hazard of the specified 
post-release/discharge drug use in the VTC group 
versus the CDDC group. 
The unadjusted hazard ratios are estimated 
using Cox regression models where group is the 
only covariate and no additional weighting is used.
The adjusted Cox models additionally adjust for
measurements of opiate cravings 
and ambivalence, recognition, 
and taking steps towards change and apply 
inverse propensity of treatment weights. 
The adjusted, time varying group effects modify the adjusted model by 
incorporating a time-varying (non-proportional) hazard for the group, 
implemented as an interaction between group and
the logarithm of time beyond 50 days post-release/discharge.
The row label reads CDDC imputation model + VTC imputation model. 
Each table uses a separate event definition as described in the subtitle.^2n;

ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}&text";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}VTC: Voluntary Treatment Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}CDDC: Compulsory Drug Detention Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SR: Stochastic regression imputation";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV90: Single value imputation of length of inpatient stay (90 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV30: Single value imputation of length of inpatient stay (30 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV14: Single value imputation of time of first outcome measurement (14 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV0: Single value imputation of time of first outcome measurement (0 days)";

%end;

%MEND;

****************************************************************************;

%let data = cox1;
%let site = site3;
%let ref =LAST;
%let type_list = A B C H;

%let var = opi;
%master_HR_estimates(&var,&type_list,&data,&site,&ref);

%let var = amp;
%master_HR_estimates(&var,&type_list,&data,&site,&ref);

%let var = any;
%master_HR_estimates(&var,&type_list,&data,&site,&ref);

%let type_list = A B C G;
%let var = mtd;
%master_HR_estimates(&var,&type_list,&data,&site,&ref);

options nodate nonumber orientation=landscape;
%let file = "&out_dir.\CDDC Urine Supp Cox AHR &var &date..rtf";
ods rtf file = &file Style = custom_narrow startpage = Yes BODYTITLE;

%let type_list = A B C H;
%let var = opi;
%let title1 = Table S23$#. @ hazard ratios for opiate use.;
%print_HR_estimates(&var,&type_list,&title1);

%let var = amp;
%let title1 = Table S24$#. @ hazard ratios for amphetamine use.;
%print_HR_estimates(&var,&type_list,&title1);

%let var = any;
%let title1 = Table S25$#. @ hazard ratios for illicit-drug use.;
%print_HR_estimates(&var,&type_list,&title1);

%let var = mtd;
%let title1 = Table S26$#. @ hazard ratios for linkage to methadone.;
%let type_list = A B C;
%print_HR_estimates(&var,&type_list,&title1);

%let var = mtd;
%let title1 = Table S26d#. @ hazard ratios for methadone retention.;
%let type_list = G;
%print_HR_estimates(&var,&type_list,&title1);


%let text = 
Note: For each stage 1 imputation combination, 
hazard ratios and 95% confidence intervals 
comparing the hazard of the specified 
post-release/discharge methadone linkage or retention in the VTC group 
versus the CDDC group. 
The unadjusted hazard ratios are estimated 
using Cox regression models where group is the 
only covariate and no additional weighting is used.
The adjusted Cox models additionally adjust for
measurements of opiate cravings 
and ambivalence, recognition, 
and taking steps towards change and apply 
inverse propensity of treatment weights. 
The adjusted, time varying group effects modify the adjusted model by 
incorporating a time-varying (non-proportional) hazard for the group, 
implemented as an interaction between group and
the logarithim of time beyond 50 days post-release/discharge.
The row label reads CDDC imputation model + VTC imputation model. 
Each table uses a separate event definition as described in the subtitle.^2n;

ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}&text";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}VTC: Voluntary Treatment Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}CDDC: Compulsory Drug Detention Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SR: Stochastic regression imputation";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV90: Single value imputation of length of inpatient stay (90 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV30: Single value imputation of length of inpatient stay (30 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV14: Single value imputation of time of first outcome measurement (14 days)";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}SV0: Single value imputation of time of first outcome measurement (0 days)";


ods rtf close;

****************************************************************************;
****************************************************************************;
*****************                 Table 3                     *************;
****************************************************************************;
****************************************************************************;

proc format;

value group_ 	1 = "Unadjusted"
				2 = "Adjusted"
				3 = "Adjusted, time-varying group effect";

picture lc (round) /* Lower limit requires a left parenthesis and a comma on the right */
 			low-high = '009.99,'
 			( prefix = '(' );
picture uc (round)
 			low-high = '009.99)'; /* Upper limit requires only a right parenthesis */

value $contrast_ "1" = "Group"
				"2" = "Group (0 days < T < 50 days)"
				"3" = "Group (T = 90 days)"
				"4" = "Group (T = 180 days)"
				"5" = "Group (T = 270 days)"
				"6" = "Group (T = 365 days)"
				"7" = "Opiate cravings"
				"8" = "Recognition of change"
				"9" = "Ambivalence towards change"
				"99" = "Taking steps towards change";
run;

data loc.cox2;
	set loc.cox1;
	if combination = 1;
run;

%let var = opi;
%let type = A;
%let site = site3;
%let ref = LAST;
%let data = cox2;

%let T = T&type._&var;
%let Y = Y&type._&var;

ods exclude all;

proc phreg data = loc.&data;
	by combination _IMPUTATION_;
	class &site (ref=&ref);
	format &site site3b_.;
	model &T*&Y(0)= &site /ties=efron;
	ods output  ContrastEstimate=est1;
	contrast "1" &site 1 /ESTIMATE=Parm;
run;

proc phreg data = loc.&data;
	by combination _IMPUTATION_;
	class &site (ref=&ref);
	format &site site3b_.;
	weight wt2;
	model &T*&Y(0)= &site &var3 /ties=efron;
	ods output  ContrastEstimate=est2;
	contrast "1" &site 1 /ESTIMATE=Parm;
	contrast "7" cs01_craving 1 /ESTIMATE=Parm;
	contrast "8" so_rec_d 10 /ESTIMATE=Parm;
	contrast "9" so_amb_d 10 /ESTIMATE=Parm;
	contrast "99" so_ts_d 10 /ESTIMATE=Parm;
run;

proc phreg data = loc.&data;
	by combination _IMPUTATION_;
	class &site (ref=&ref);	
	format &site site3b_.;
	weight wt2;
	model &T*&Y(0)= &site site50lt &var3 /ties=efron;
	if &T > 50 then site50lt = &site*log((&T-50));
	else site50lt = 0;
	ods output  ContrastEstimate=est3;
	contrast "2" &site 1 /ESTIMATE=Parm;
	contrast "3" &site 1 site50lt 3.688 /ESTIMATE=Parm;
	contrast "4" &site 1 site50lt 4.868 /ESTIMATE=Parm;
	contrast "5" &site 1 site50lt 5.394 /ESTIMATE=Parm;
	contrast "6" &site 1 site50lt 5.753 /ESTIMATE=Parm;
	contrast "7" cs01_craving 1 /ESTIMATE=Parm;
	contrast "8" so_rec_d 10 /ESTIMATE=Parm;
	contrast "9" so_amb_d 10 /ESTIMATE=Parm;
	contrast "99" so_ts_d 10 /ESTIMATE=Parm;
run;

proc sort data=est1;
	by contrast;
proc sort data=est2;
	by contrast;
proc sort data=est3;
	by contrast;
data est_all;
	length contrast $30.;
	set est1 (in = in1) est2 (in = in2) est3 (in = in3);
	by contrast;
	if in1 then group = 1;
	if in2 then group = 2;
	if in3 then group = 3;
run;

proc sort data=est_all;
	by group contrast _imputation_;
run;

proc mianalyze data=est_all;
	title;
	by group contrast;
	modeleffects estimate;
	stderr stderr;
	ods output ParameterEstimates=parms;
run;

ods exclude none;

data parms2;
	set parms;

	if Parm = "estimate";

	HR = exp(-Estimate);
	LCL = Estimate - 1.96*StdErr;
	UCL = Estimate + 1.96*StdErr;
	lowerCI = exp(-UCL);
	upperCI = exp(-LCL);
	CI = "(" || put(lowerCI,z5.3) || " - " || put(upperCI,z5.3) || ")";

	keep contrast group HR lowerCI upperCI;
run;

options nodate nonumber orientation=landscape;
%let file = "&out_dir./CDDC Urine Table 3 &date..rtf";
ods rtf file = &file Style = custom_narrow startpage = Yes BODYTITLE;

title "Table 3. Hazard Ratios for Study Variables.";
proc tabulate data= parms2;
	class contrast group;
	var HR lowerCI upperCI;
	table 	Contrast,
		group=' ' *(HR*max='' * [style=[just=C cellwidth=75] f=z5.3]
 					lowerCI*max="" * [style=[just=R cellwidth=65 ] f=lc.]
 					upperCI*max="" * [style=[just=L cellwidth=65 ] f=uc.]);	
	format group group_. contrast $contrast_.;
	label HR="Hazard Ratio";
	label contrast="Variable";
	label lowerCI ='95%';
	label upperCI ="C.I.";
run;
title;

%let text = 
Note: Table provides the hazard ratios and 95% confidence intervals 
comparing the hazard of the specified 
post-release/discharge opirate relapse in the VTC group 
versus the CDDC group. 
The unadjusted hazard ratios are estimated 
using Cox regression models where group is the 
only covariate and no additional weighting is used.
The adjusted Cox models additionally adjust for
measurements of opiate cravings 
and ambivalence, recognition, 
and taking steps towards change and apply 
inverse propensity of treatment weights. 
The adjusted, time varying group effects modify the adjusted model by 
incorporating a time-varying (non-proportional) hazard for the group, 
implemented as an interaction between group and
the logarithm of time beyond 50 days post-release/discharge. 
Stochastic regression for missing discharge and inpatient times are used
missing outcome measurements are ignored. ^2n;

ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}&text";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}VTC: Voluntary Treatment Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}CDDC: Compulsory Drug Detention Center";
ods text = "^S={LEFTMARGIN = 1in RIGHTMARGIN = 1in}C.I.: Confidence Interval";

ods rtf close;


****************************************************************************;
****************************************************************************;
*********************   Remove old datasets               ******************;
****************************************************************************;
****************************************************************************;


proc datasets lib=loc nolist;       
	delete cox1 temp3;   
run;
quit;

****************************************************************************;
****************************************************************************;
*****************                   End                        *************;
****************************************************************************;
****************************************************************************;

%MEND cox_analysis;

